unit Functions;

{}

interface


	uses
		QuickDraw, Palettes, PrintTraps, globals, Utilities, Graphics, File1, File2, Analysis, Camera, Lut;


	procedure ApplyTable (var table: LookupTable);
	procedure ApplyLookupTable;
	procedure MakeBinary;
	procedure Filter (ftype: FilterType; pass: integer; var table: FateTable);
	procedure PhotoMode;
	function AllSameSize: boolean;
	procedure EnhanceContrast;
	procedure EqualizeHistogram;
	procedure Convolve (name: str255; RefNum: integer);
	procedure ConvolveUsingText;
	procedure PlotSurface;
	procedure MakeSkeleton;
	procedure DoErosion;
	procedure DoDilation;
	procedure DoOpening;
	procedure DoClosing;
	procedure SetBinaryCount;
	procedure SetIterations;
	procedure ChangeValues (v1, v2, v3: integer);
	procedure DoPropagate (MenuItem: integer);
	procedure DoArithmetic (MenuItem: integer; constant: extended);
	procedure NewPlotSurface;
	procedure AutoThreshold;
	procedure AutoDensitySlice;
	procedure FixColors;
	procedure DoImageMath;


implementation

	const
		MaxW = 4000;
		Src1Item = 7;
		Src2Item = 8;
		OpItem = 9;

	type
		ktype = array[0..MaxW] of integer;
		SortArray = array[1..9] of integer;

	var
		PixelsRemoved: LongInt;
		Src1PicNum, Src2PicNum: integer;

	procedure ApplyTableToLine (data: ptr; var table: LookupTable; width: LongInt);
{$IFC false}
		var
			line: LinePtr;
			i: integer;
	begin
		line := LinePtr(data);
		for i := 0 to width - 1 do
			Line^[i] := table[Line^[i]];
	end;
{$ENDC}

{a0 = data}
{a1 = lookup table}
{d0 = width }
{d1 = pixel value}
inline
	$4E56, $0000, {  link a6,#0}
	$48E7, $C0C0, {  movem.l a0-a1/d0-d1,-(sp)}
	$206E, $000C, {  move.l 12(a6),a0}
	$226E, $0008, {  move.l 8(a6),a1}
	$202E, $0004, {  move.l 4(a6),d0}
	$5380,       {  subq.l #1,d0}
	$4281,       {  clr.l d1}
	$1210,       {L move.b (a0),d1}
	$10F1, $1000, {  move.b 0(a1,d1.w),(a0)+}
	$51C8, $FFF8, {  dbra d0,L}
	$4CDF, $0303, {  movem.l (sp)+,a0-a1/d0-d1}
	$4E5E,       {  unlk a6}
	$DEFC, $000C; {  add.w #12,sp}


procedure PutLineUsingMask (h, v, count: integer; var line: LineType);
	var
		aLine, MaskLine: LineType;
		i: integer;
		SaveInfo: InfoPtr;
begin
	if count > MaxLine then
		count := MaxLine;
	GetLine(h, v, count, aline);
	SaveInfo := Info;
	Info := UndoInfo;
	GetLine(h, v, count, MaskLine);
	for i := 0 to count - 1 do
		if MaskLine[i] = BlackIndex then
			aLine[i] := line[i];
	info := SaveInfo;
	PutLine(h, v, count, aLine);
end;


procedure ApplyTable; {(var table: LookupTable)}
	var
		width, NumberOfLines, i, hloc, vloc: integer;
		offset: LongInt;
		p: ptr;
		UseMask: boolean;
		TempLine: LineType;
		AutoSelectAll: boolean;
begin
	if NotInBounds then
		exit(ApplyTable);
	StopDigitizing;
	AutoSelectAll := not Info^.RoiShowing;
	if AutoSelectAll then
		SelectAll(false);
	if TooWide then
		exit(ApplyTable);
	ShowWatch;
	with info^.RoiRect, info^ do begin
			if RoiType <> RectRoi then
				UseMask := SetupMask
			else
				UseMask := false;
			SetupUndoFromClip;
			WhatToUndo := UndoTransform;
			offset := LongInt(top) * BytesPerRow + left;
			if UseMask then
				p := @TempLine
			else
				p := ptr(ord4(PicBaseAddr) + offset);
			width := right - left;
			NumberOfLines := bottom - top;
			hloc := left;
			vloc := top;
		end;
	if width > 0 then
		for i := 1 to NumberOfLines do
			if UseMask then begin
					GetLine(hloc, vloc, width, TempLine);
					ApplyTableToLine(p, table, width);
					PutLineUsingMask(hloc, vloc, width, TempLine);
					vloc := vloc + 1
				end
			else begin
					ApplyTableToLine(p, table, width);
					p := ptr(ord4(p) + info^.BytesPerRow);
				end;
	with info^ do begin
			UpdateScreen(RoiRect);
			Info^.changes := true;
		end;
	SetupRoiRect;
	if AutoSelectAll then
		KillRoi;
end;


function DoApplyTableDialogBox: boolean;
	const
		Button1 = 3;
		Button2 = 4;
		Button3 = 5;
		Button4 = 6;
	var
		mylog: DialogPtr;
		item: integer;
		SaveA, SaveB: boolean;

	procedure SetButtons;
	begin
		SetDialogItem(mylog, Button1, ord(ThresholdToForeground));
		SetDialogItem(mylog, Button2, ord(not ThresholdToForeground));
		SetDialogItem(mylog, Button3, ord(NonThresholdToBackground));
		SetDialogItem(mylog, Button4, ord(not NonThresholdToBackground));
	end;

begin
	InitCursor;
	SaveA := ThresholdToForeground;
	SaveB := NonThresholdToBackground;
	mylog := GetNewDialog(40, nil, pointer(-1));
	SetButtons;
	OutlineButton(MyLog, ok, 16);
	repeat
		ModalDialog(nil, item);
		if (item = Button1) or (item = button2) then begin
				ThresholdToForeground := not ThresholdToForeground;
				SetButtons;
			end;
		if (item = Button3) or (item = button4) then begin
				NonThresholdToBackground := not NonThresholdToBackground;
				SetButtons;
			end;
	until (item = ok) or (item = cancel);
	DisposDialog(mylog);
	if item = cancel then begin
			ThresholdToForeground := SaveA;
			NonThresholdToBackground := SaveB;
			DoApplyTableDialogBox := false
		end
	else
		DoApplyTableDialogBox := true;
end;


procedure ApplyLookupTable;
	var
		table: LookupTable;
		ConvertingColorPic, GrayScaleImage: boolean;
		i: integer;
begin
	with info^ do begin
			GrayScaleImage := (LUTMode = Grayscale) or (LUTMode = CustomGrayscale);
			ConvertingColorPic := not GrayScaleImage and not DensitySlicing;
			if ConvertingColorPic then
				KillRoi;
			if DensitySlicing and (not macro) then begin
					if not DoApplyTableDialogBox then
						exit(ApplyLookupTable);
				end;
			if thresholding then
				BinaryPic := true;
			GetLookupTable(table);
			if GrayscaleImage or ConvertingColorPic then
				ResetGrayMap;
			ApplyTable(table);
			if ConvertingColorPic then
				WhatToUndo := NothingToUndo;
			if DensityCalibrated then begin
					DensityCalibrated := false;
					for i := 0 to 255 do
						cvalue[i] := i;
				end;
		end; {with}
end;


procedure MakeBinary;
	var
		table: LookupTable;
		SaveBackground, SaveForeground, i: integer;
begin
	with info^ do begin
			if DensitySlicing then begin
					ThresholdToForeground := true;
					NonThresholdToBackground := true;
					SaveBackground := BackgroundIndex;
					SaveForeground := ForegroundIndex;
					BackgroundIndex := WhiteIndex;
					ForegroundIndex := BlackIndex;
					GetLookupTable(table);
					ResetGrayMap;
					ApplyTable(table);
					BackgroundIndex := SaveBackground;
					ForegroundIndex := SaveForeground;
					BinaryPic := true;
				end
			else if Thresholding then begin
					for i := 0 to 255 do
						if i < ColorStart then
							table[i] := WhiteIndex
						else
							table[i] := BlackIndex;
					ResetGrayMap;
					ApplyTable(table);
					BinaryPic := true;
				end
			else
				PutMessage('Sorry, but you must be thresholding or density slicing to use Make Binary.');
		end;
end;


{$IFC false}
function FindMedian (var a: SortArray): integer;
  {Finds the 5th largest of 9 values}
	var
		i, j, mj, max: integer;
begin
	for i := 1 to 4 do begin
			max := 0;
			mj := 1;
			for j := 1 to 9 do
				if a[j] > max then begin
						max := a[j];
						mj := j;
					end;
			a[mj] := 0;
		end;
	max := 0;
	for j := 1 to 9 do
		if a[j] > max then
			max := a[j];
	FindMedian := max;
end;
{$ENDC}

function FindMedian (var a: sortArray): integer;
   {In-line code contributed by Edward J. Huff(huff@mcclbo.med.nyu.edu).}
   {Assember source with comments and a test program are available by anonymous}
   {ftp from zippy.nimh.nih.gov, in the /pub/nih-image/documents directory.}
inline
	$205F, $48E7, $1F00, $4C98, $00FF, $B041, $6502, $C340,{}
	$B443, $6502, $C742, $B243, $6504, $C540, $C741, $B845,{}
	$6502, $CB44, $BC47, $6502, $CF46, $BA47, $6504, $CD44,{}
	$CF45, $B245, $6508, $CF43, $CD42, $CB41, $C940, $3E10,{}
	$BC47, $6502, $CF46, $BA47, $6504, $CD44, $CF45, $B245,{}
	$6508, $CF43, $CD42, $CB41, $C940, $B246, $6534, $B242,{}
	$6514, $B244, $6504, $3001, $6062, $B644, $6504, $3004,{}
	$605A, $3003, $6056, $B444, $650C, $B445, $6504, $3005,{}
	$604A, $3002, $6046, $B644, $6504, $3004, $603E, $3003,{}
	$603A, $B645, $6504, $C942, $CB43, $B846, $651C, $B644,{}
	$650C, $B444, $6504, $3002, $6022, $3004, $601E, $B646,{}
	$6504, $3003, $6016, $3006, $6012, $B646, $6508, $B446,{}
	$65F4, $3002, $6006, $B644, $65E0, $3003, $4CDF, $00F8,{}
	$3E80;


procedure Filter (ftype: FilterType; pass: integer; var table: FateTable);
	const
		PixelsPerUpdate = 5000;
	var
		row, width, r1, r2, r3, c, value, error, sum, center: integer;
		tmp, mark, NewMark, LinesPerUpdate, LineCount: integer;
		t1, t2, t3, t4: integer;
		MaskRect, frame, trect: rect;
		WhitePixel1: integer;
		L1: LineType;
		WhitePixel2: integer;
		L2: LineType;
		WhitePixel3: integer;
		L3, result: LineType;
		pt: point;
		a: SortArray;
		AutoSelectAll, UseMask, BinaryFilter: boolean;
		L, T, R, B, index, code, FirstRow, LastRow: integer;
		StartTicks: LongInt;
begin
	if NotinBounds then
		exit(Filter);
	StopDigitizing;
	AutoSelectAll := not Info^.RoiShowing;
	if AutoSelectAll then
		with info^ do begin
				SelectAll(false);
				SetPort(wptr);
				PenNormal;
				PenPat(pat[PatIndex]);
				FrameRect(wrect);
			end;
	if TooWide then
		exit(Filter);
	ShowWatch;
	if info^.RoiType <> RectRoi then
		UseMask := SetupMask
	else
		UseMask := false;
	if pass = 0 then begin
			SetupUndoFromClip;
			ShowMessage(CmdPeriodToStop);
			WhatToUndo := UndoFilter;
		end;
	frame := info^.RoiRect;
	StartTicks := TickCount;
	BinaryFilter := ftype in [Erosion, Dilation, OutlineFilter, Skeletonize];
	with frame, Info^ do begin
			changes := true;
			RoiShowing := false;
			width := right - left;
			LinesPerUpdate := PixelsPerUpdate div width;
			if ftype = ReduceNoise then
				LinesPerUpdate := LinesPerUpdate div 3;
			if BinaryFilter then begin
					FirstRow := top;
					LastRow := bottom - 1;
					WhitePixel1 := WhiteIndex;
					WhitePixel2 := WhiteIndex;
					WhitePixel3 := WhiteIndex;
					if width < MaxLine then begin
							L1[width] := WhiteIndex;
							L2[width] := WhiteIndex;
							L3[width] := WhiteIndex;
						end;
				end
			else begin
					FirstRow := top + 1;
					LastRow := bottom - 2;
				end;
			GetLine(left, FirstRow - 1, width, L2);
			GetLine(left, FirstRow, width, L3);
			Mark := RoiRect.top;
			LineCount := 0;
			for row := FirstRow to LastRow do begin
       {Move Convolution Window Down}
					BlockMove(@L2, @L1, width);
					BlockMove(@L3, @L2, width);
					GetLine(left, row + 1, width, L3);
       {Process One Row}
					case ftype of
						EdgeDetect: 
							for c := 1 to width - 2 do begin
									t1 := L1[c - 1] + L1[c] + L1[c + 1] - L3[c - 1] - L3[c] - L3[c + 1];
									t1 := abs(t1);
									t2 := L1[c + 1] + L2[c + 1] + L3[c + 1] - L1[c - 1] - L2[c - 1] - L3[c - 1];
									t2 := abs(t2);
									if t1 > t2 then
										tmp := t1
									else
										tmp := t2;
									if OptionKeyWasDown then begin
											if tmp > 255 then
												tmp := 255;
											if tmp < 0 then
												tmp := 0;
										end
									else if tmp > 35 then
										tmp := 255
									else
										tmp := 0;
									result[c] := tmp;
								end;
						ReduceNoise:  {Median Filter}
							for c := 1 to width - 2 do begin
									a[1] := L1[c - 1];
									a[2] := L1[c];
									a[3] := L1[c + 1];
									a[4] := L2[c - 1];
									a[5] := L2[c];
									a[6] := L2[c + 1];
									a[7] := L3[c - 1];
									a[8] := L3[c];
									a[9] := L3[c + 1];
									result[c] := FindMedian(a);
								end;
						Dither:  {Floyd-Steinberg Algorithm}
							for c := 1 to width - 2 do begin
									value := L2[c];
									if value < 128 then begin
											result[c] := 0;
											error := -value;
										end
									else begin
											result[c] := 255;
											error := 255 - value
										end;
									tmp := L2[c + 1];              {A}
									tmp := tmp - (7 * error) div 16;
									if tmp < 0 then
										tmp := 0;
									if tmp > 255 then
										tmp := 255;
									L2[c + 1] := tmp;
									tmp := L3[c + 1];              {B}
									tmp := tmp - error div 16;
									if tmp < 0 then
										tmp := 0;
									if tmp > 255 then
										tmp := 255;
									L3[c + 1] := tmp;
									tmp := L3[c];              {C}
									tmp := tmp - (5 * error) div 16;
									if tmp < 0 then
										tmp := 0;
									if tmp > 255 then
										tmp := 255;
									L3[c] := tmp;
									tmp := L3[C - 1];                {D}
									tmp := tmp - (3 * error) div 16;
									if tmp < 0 then
										tmp := 0;
									if tmp > 255 then
										tmp := 255;
									L3[C - 1] := tmp;
								end;
						UnweightedAvg: 
							for c := 1 to width - 2 do begin
									tmp := (L1[C - 1] + L1[c] + L1[c + 1] + L2[c - 1] + L2[c] + L2[c + 1] + L3[c - 1] + L3[c] + L3[c + 1]) div 9;
									if tmp > 255 then
										tmp := 255;
									if tmp < 0 then
										tmp := 0;
									result[c] := tmp;
								end;
						WeightedAvg: 
							for c := 1 to width - 2 do begin
									tmp := (L1[c - 1] + L1[c] + L1[c + 1] + L2[c - 1] + L2[c] * 4 + L2[c + 1] + L3[c - 1] + L3[c] + L3[c + 1]) div 12;
									if tmp > 255 then
										tmp := 255;
									if tmp < 0 then
										tmp := 0;
									result[c] := tmp;
								end;
						fsharpen: 
							for c := 1 to width - 2 do begin
									if OptionKeyWasDown then
										tmp := L2[c] * 9 - L1[c - 1] - L1[c] - L1[c + 1] - L2[c - 1] - L2[c + 1] - L3[c - 1] - L3[c] - L3[c + 1]
									else begin
											tmp := L2[c] * 12 - L1[c - 1] - L1[c] - L1[c + 1] - L2[c - 1] - L2[c + 1] - L3[c - 1] - L3[c] - L3[c + 1];
											tmp := tmp div 4;
										end;
									if tmp > 255 then
										tmp := 255;
									if tmp < 0 then
										tmp := 0;
									result[c] := tmp;
								end;
						fshadow: 
							for c := 1 to width - 2 do begin
									tmp := L2[c + 1] + L2[c + 1] + L3[c] + L3[c + 1] * 2 - L1[c - 1] * 2 - L1[c] - L2[c - 1];
									if tmp > 255 then
										tmp := 255;
									if tmp < 0 then
										tmp := 0;
									result[c] := tmp;
								end;
						Erosion: 
							for c := 0 to width - 1 do begin
									center := L2[c];
									if center = BlackIndex then begin
											sum := L1[c - 1] + L1[c] + L1[c + 1] + L2[c - 1] + L2[c + 1] + L3[c - 1] + L3[c] + L3[c + 1];
											if (2040 - sum) >= BinaryThreshold then
												center := WhiteIndex;
										end;
									result[c] := center;
								end;
						Dilation: 
							for c := 0 to width - 1 do begin
									center := L2[c];
									if center = WhiteIndex then begin
											sum := L1[c - 1] + L1[c] + L1[c + 1] + L2[c - 1] + L2[c + 1] + L3[c - 1] + L3[c] + L3[c + 1];
											if sum >= BinaryThreshold then
												center := BlackIndex;
										end;
									result[c] := center;
								end;
						OutlineFilter: 
							for c := 0 to width - 1 do begin
									center := L2[c];
									if center = BlackIndex then begin
											if (L2[c - 1] = WhiteIndex) or (L1[c] = WhiteIndex) or (L2[c + 1] = WhiteIndex) or (L3[c] = WhiteIndex) then
												center := BlackIndex
											else
												center := WhiteIndex;
										end;
									result[c] := center;
								end;

						Skeletonize: 
							for c := 0 to width - 1 do begin
									center := L2[c];
									if center = BlackIndex then begin
											index := 0;
											if L1[c - 1] = BlackIndex then
												index := bor(index, 1);
											if L1[c] = BlackIndex then
												index := bor(index, 2);
											if L1[c + 1] = BlackIndex then
												index := bor(index, 4);
											if L2[c + 1] = BlackIndex then
												index := bor(index, 8);
											if L3[c + 1] = BlackIndex then
												index := bor(index, 16);
											if L3[c] = BlackIndex then
												index := bor(index, 32);
											if L3[c - 1] = BlackIndex then
												index := bor(index, 64);
											if L2[c - 1] = BlackIndex then
												index := bor(index, 128);
											code := table[index];
											if odd(pass) then begin
													if (code = 2) or (code = 3) then begin
															center := WhiteIndex;
															PixelsRemoved := PixelsRemoved + 1;
														end;
												end
											else begin {even pass}
													if (code = 1) or (code = 3) then begin
															center := WhiteIndex;
															PixelsRemoved := PixelsRemoved + 1;
														end;
												end;
										end; {if}
									result[c] := center;
								end; {for}
					end; {case}
					if not BinaryFilter then begin
							result[0] := L2[0];
							result[width - 1] := L2[width - 1];
						end;
					if UseMask then
						PutLineUsingMask(left, row, width, result)
					else
						PutLine(left, row, width, result);
					LineCount := LineCount + 1;
					if LineCount = LinesPerUpdate then begin
							pt.h := RoiRect.left;
							pt.v := row + 1;
							NewMark := pt.v;
							with RoiRect do
								SetRect(MaskRect, left, mark, right, NewMark);
							UpdateScreen(MaskRect);
							LineCount := 0;
							Mark := NewMark;
							if magnification > 1.0 then
								Mark := Mark - 1;
							if CommandPeriod then begin
									UpdatePicWindow;
									beep;
									PixelsRemoved := 0;
									if AutoSelectAll then
										KillRoi;
									exit(filter)
								end;
						end;
				end; {for row:=...}
			trect := frame;
			InsetRect(trect, 1, 1);
			ShowTime(StartTicks, trect, '');
		end; {with}
	if LineCount > 0 then begin
			with frame do
				SetRect(MaskRect, left, mark, right, bottom);
			UpdateScreen(MaskRect)
		end;
	SetupRoiRect;
	if AutoSelectAll then
		KillRoi;
end;


procedure PhotoMode;
{Erases the screen to the background color and then redraws}
{the contents of the active image window . }
	var
		tPort: GrafPtr;
		event: EventRecord;
		WinRect: rect;
		SaveVisRgn: rgnHandle;
begin
	with info^ do begin
			KillRoi;
			if OptionKeyWasDown then begin {Move window up to top of screen.}
					GetWindowRect(wptr, WinRect);
					MoveWindow(wptr, WinRect.left, 0, false);
				end;
			with wptr^ do begin
					SaveVisRgn := visRgn;
					visRgn := NewRgn;
					RectRgn(visRgn, ScreenBits.Bounds);
				end;
			FlushEvents(EveryEvent, 0);
			GetPort(tPort);
			EraseScreen;
			UpdatePicWindow;
			repeat
			until WaitNextEvent(mDownMask + KeyDownMask, Event, 5, nil);
			with wptr^ do begin
					DisposeRgn(visRgn);
					visRgn := SaveVisRgn;
				end;
			RestoreScreen;
			SetPort(tPort);
			FlushEvents(EveryEvent, 0);
			if OptionKeyWasDown then begin
					MoveWindow(wptr, WinRect.left, WinRect.top, false);
				end;
		end;
end;


function AllSameSize: boolean;
{Returns true if all currently open Images have the same dimensions.}
	var
		i: integer;
		SameSize: Boolean;
		TempInfo: InfoPtr;
begin
	if nPics = 0 then begin
			AllSameSize := false;
			exit(AllSameSize);
		end;
	SameSize := true;
	for i := 1 to nPics do begin
			TempInfo := pointer(WindowPeek(PicWindow[i])^.RefCon);
			SameSize := SameSize and EqualRect(Info^.PicRect, TempInfo^.PicRect);
		end;
	AllSameSize := SameSize;
end;


procedure EnhanceContrast;
	var
		AutoSelectAll: boolean;
		min, max, i, threshold: integer;
		found, SaveRedirectFlag: boolean;
		sum: LongInt;
begin
	with info^ do
		if LUTMode = ColorLUT then begin
				PutMessage('Sorry, but you can not contrast enhance true color images.');
				exit(EnhanceContrast)
			end;
	if NotInBounds or (ClipBuf = nil) then
		exit(EnhanceContrast);
	StopDigitizing;
	AutoSelectAll := not Info^.RoiShowing;
	if AutoSelectAll then
		SelectAll(false);
	SaveRedirectFlag := RedirectSampling;
	RedirectSampling := false;
	if info^.RoiType = RectRoi then
		GetRectHistogram
	else
		GetHistogram;
	RedirectSampling := SaveRedirectFlag;
	sum := 0;
	for i := 0 to 255 do
		sum := sum + histogram[i];
	threshold := sum div 5000;
	i := -1;
	repeat
		i := i + 1;
		found := histogram[i] > threshold;
	until found or (i = 255);
	min := i;
	i := 256;
	repeat
		i := i - 1;
		found := histogram[i] > threshold;
	until found or (i = 0);
	max := i;
	if max > min then
		with info^ do begin
				SetupLutUndo;
				if isGrayScaleLUT then
					LUTMode := grayscale;
				ColorStart := min;
				ColorEnd := max;
				DrawMap;
				UpdateLUT;
				changes := true;
				IdentityFunction := false;
			end;
	if AutoSelectAll then
		KillRoi;
end;


procedure EqualizeHistogram;
	var
		AutoSelectAll, SaveRedirectFlag: boolean;
		i, sum, v: integer;
		isum: LongInt;
		ScaleFactor: extended;
begin
	with info^ do
		if (LUTMode <> GrayScale) and (LutMode <> CustomGrayscale) then begin
				PutMessage('Sorry, but you can only do histogram equalization on grayscale images.');
				exit(EqualizeHistogram)
			end;
	if NotInBounds or (ClipBuf = nil) then
		exit(EqualizeHistogram);
	StopDigitizing;
	AutoSelectAll := not Info^.RoiShowing;
	if AutoSelectAll then
		SelectAll(false);
	SaveRedirectFlag := RedirectSampling;
	RedirectSampling := false;
	if info^.RoiType = RectRoi then
		GetRectHistogram
	else
		GetHistogram;
	RedirectSampling := SaveRedirectFlag;
	FindThresholdingMode;
	ComputeResults;
	isum := 0;
	for i := 0 to 255 do
		isum := isum + histogram[i];
	ScaleFactor := 255.0 / isum;
	sum := 0;
	with info^ do begin
			SetupLutUndo;
			for i := 255 downto 0 do
				with cTable[i].rgb do begin
						sum := round(sum + histogram[i] * ScaleFactor);
						if sum > 255 then
							sum := 255;
						v := sum * 256;
						red := v;
						green := v;
						blue := v;
					end;
			LoadLUT(cTable);
			LUTMode := CustomGrayscale;
			SetupPseudocolor;
			changes := true;
			DrawMap;
			IdentityFunction := false;
		end; {with info}
	if AutoSelectAll then
		KillRoi;
end;


procedure GetKernel (var kernel: ktype; var n: integer; var name: str255; RefNum: integer);
	var
		rLine: rLineType;
		i, count, nValues, nRows: integer;
begin
	count := 0;
	nRows := 0;
	InitTextInput(name, RefNum);
	while not TextEof and (nRows <= 63) do begin
			GetLineFromText(rLine, nValues);
			if count <> 0 then
				nRows := nRows + 1;
			if nRows = 1 then
				n := nValues;
			for i := 1 to nValues do begin
					count := count + 1;
					kernel[count - 1] := round(rLine[i]);
				end;
		end;
	if count <> (n * n) then
		n := 0;
end;


procedure DoOnePixel (nLess1, BytesPerLine: integer; corner: LongInt; var sum: LongInt; var kernel: ktype);
{$IFC false}
	var
		row, column, k: integer;
		pp: ptr;
begin
	k := 0;
	sum := 0;
	for row := 0 to nless1 do begin
			corner := corner + BytesPerLine;
			pp := ptr(corner);
			for column := 0 to nless1 do begin
					sum := sum + band(pp^, 255) * kernel[k];
					k := k + 1;
					pp := ptr(ord(pp) + 1);
				end;
		end;
end;
{$ENDC}

{a0=^corner/^sum}
{a1=^kernel}
{a2=^pixels}

{d0=n-1}
{d1=BytesPerLine}
{d2=sum}
{d3=n-1(outer loop)}
{d4=n-1(inner loop)}
{d5=temp}

inline
	$4E56, $0000, {  link	a6,#0}
	$48E7, $FCE0,  {  movem.l	a0-a2/d0-d5,-(sp)}
	$4280,              {  clr.l	d0}
	$302E, $0012, {  move.w	18(a6),d0}
	$4281,              {  clr.l	d1}
	$322E, $0010, {  move.w	16(a6),d1}
	$206E, $000C, {  movea.l	12(a6),a0}
	$226E, $0004, {  movea.l	4(a6),a1}

	$4282,             {  clr.l	d2}
	$2600,             {  move.l	d0,d3}

	$D1C1,             {A adda.l	d1,a0}
	$2448,            {  move.l	a0,a2}
	$2800,            {  move.l	d0,d4}
	$4285,            {B clr.l	d5                   (2)}
	$1A1A,             {  move.b	(a2)+,d5    (6) }
	$CBD9,             {  muls	(a1)+,d5     (29!)}
	$D485,             {  add.l	d5,d2          (2)}
	$51CC, $FFF6, {  dbra	d4,B                (6)}
	$51CB, $FFEC, {  dbra	d3,A}

	$206E, $0008, {  move.l	8(a6),a0}
	$2082,              {  move.l	d2,(a0)}
	$4CDF, $073F, {  movem.l	(sp)+,a0-a2/d0-d5}
	$4E5E,              {  unlk	a6}
	$DEFC, $0010; {  add.w	#16,sp}



procedure DoConvolution (var kernel: ktype; n: integer);
	const
		skip = 7;
	var
		row, width, column, value, error: integer;
		margin, i, nless1: integer;
		frame, MaskRect, tRect: rect;
		AutoSelectAll, ScalingNeeded: boolean;
		SrcCenter, DstCenter, sum, max, offset, wsum, cscale, StartTicks: LongInt;
		MinResult, MaxResult: LongInt;
		p: ptr;
		str, str2: str255;
		ScaleFactor: extended;
begin
	if NotinBounds or NotRectangular then
		exit(DoConvolution);
	StopDigitizing;
	AutoSelectAll := not Info^.RoiShowing;
	if AutoSelectAll then
		SelectAll(false);
	SetupUndoFromClip;
	WhatToUndo := UndoFilter;
	frame := info^.RoiRect;
	with frame, Info^ do begin
			if ((LutMode = GrayScale) or (LutMode = CustomGrayscale)) and (not IdentityFunction) then
				ApplyLookupTable;
			changes := true;
			margin := n div 2;
			if left < margin then
				left := left + margin;
			if right > (PicRect.right - margin) then
				right := right - margin;
			if top < margin then
				top := top + margin;
			if bottom > (PicRect.bottom - margin) then
				bottom := bottom - margin;
			SetPort(wptr);
			PenNormal;
			PenPat(pat[PatIndex]);
			tRect := frame;
			OffscreenToScreenRect(tRect);
			FrameRect(tRect);
			width := right - left;
			max := n * n - 1;
			wsum := 0;
			for i := 0 to max do
				wsum := wsum + kernel[i];
			NumToString(n, str);
			NumToString(wsum, str2);
			ValuesMessage := Concat(str, ' x ', str, ' kernel', cr, 'sum = ', str2, cr, cr, CmdPeriodToStop);
			ShowValues;
			if wsum <> 0 then
				cscale := wsum
			else
				cscale := 1;
			offset := -(n div 2) * BytesPerRow - BytesPerRow - n div 2;
			nless1 := n - 1;
			StartTicks := TickCount;
			str := '';
			if ScaleConvolutions then begin
					MinResult := MaxLongInt;
					MaxResult := -MaxLongInt;
					row := top;
					while row < bottom do begin
							SrcCenter := ord4(ClipBufInfo^.PicBaseAddr) + LongInt(row) * BytesPerRow + left;
							column := left;
							while column < (left + width) do begin
									DoOnePixel(nless1, BytesPerRow, SrcCenter + offset, sum, kernel);
									value := sum div cscale;
									if value < MinResult then
										MinResult := value;
									if value > MaxResult then
										MaxResult := value;
									SrcCenter := SrcCenter + skip;
									column := column + skip;
								end; {while column}
							row := row + skip;
						end; {while row...}
					ScalingNeeded := (MinResult < 0) or (MaxResult > 255);
					if ScalingNeeded then
						ScaleFactor := 253.0 / (MaxResult - MinResult)
					else
						ScaleFactor := 1.0;
					RealToString(ScaleFactor, 1, 4, str);
					str := concat('min=', long2str(MinResult), cr, 'max=', long2str(MaxResult), cr, 'scale factor= ', str);
					for row := top to bottom - 1 do begin
							SrcCenter := ord4(ClipBufInfo^.PicBaseAddr) + LongInt(row) * BytesPerRow + left;
							DstCenter := ord4(PicBaseAddr) + LongInt(row) * BytesPerRow + left;
							for column := left to left + width - 1 do begin
									DoOnePixel(nless1, BytesPerRow, SrcCenter + offset, sum, kernel);
									value := sum div cscale;
									if ScalingNeeded then begin
											if value < MinResult then
												value := MinResult;
											if value > MaxResult then
												value := MaxResult;
											value := round((value - MinResult) * ScaleFactor + 1);
										end;
									p := ptr(DstCenter);
									p^ := BAND(value, 255);
									SrcCenter := SrcCenter + 1;
									DstCenter := DstCenter + 1;
								end; {for column:=}
							SetRect(MaskRect, left, row, right, row + 1);
							UpdateScreen(MaskRect);
							if CommandPeriod then begin
									UpdatePicWindow;
									beep;
									exit(DoConvolution)
								end;
						end; {for row:=...}
				end  {Scale Convolutions}
			else
				for row := top to bottom - 1 do begin
						SrcCenter := ord4(ClipBufInfo^.PicBaseAddr) + LongInt(row) * BytesPerRow + left;
						DstCenter := ord4(PicBaseAddr) + LongInt(row) * BytesPerRow + left;
						for column := left to left + width - 1 do begin
								DoOnePixel(nless1, BytesPerRow, SrcCenter + offset, sum, kernel);
								value := sum div cscale;
								if value < MinResult then
									MinResult := value;
								if value > MaxResult then
									MaxResult := value;
								if value > 255 then
									value := 255;
								if value < 0 then
									value := 0;
								p := ptr(DstCenter);
								p^ := BAND(value, 255);
								SrcCenter := SrcCenter + 1;
								DstCenter := DstCenter + 1;
							end; {for column:=}
						SetRect(MaskRect, left, row, right, row + 1);
						UpdateScreen(MaskRect);
						if CommandPeriod then begin
								UpdatePicWindow;
								beep;
								exit(DoConvolution)
							end;
					end; {for row:=...}
			ShowTime(StartTicks, frame, str);
		end; {with}
	UpdatePicWindow;
	SetupRoiRect;
	if AutoSelectAll then
		KillRoi;
end;


procedure Convolve (name: str255; RefNum: integer);
	var
		kernel: ktype;
		n, count: integer;
begin
	if name = '' then begin
			RefNum := 0;
			if not GetTextFile(name, RefNum) then
				exit(convolve)
			else
				KernelsRefNum := RefNum;
		end;
	DisableDensitySlice;
	GetKernel(kernel, n, name, RefNum);
	count := n * n;
	UpdatePicWindow;
	if (n >= 3) and (n <= 63) then
		DoConvolution(kernel, n)
	else
		PutMessage('Kernels must be n x n square matrices with 3 <= n <= 63.');
end;


procedure ConvolveUsingText;
	var
		f: integer;
		err: OSErr;
		count: LongInt;
begin
	err := fsdelete('TempKernel', SystemRefNum);
	err := create('TempKernel', SystemRefNum, 'imag', 'TEXT');
	if err = NoErr then
		err := fsopen('TempKernel', SystemRefNum, f);
	if err <> NoErr then begin
			putmessage('Unable to open temporary file.');
			exit(ConvolveUsingText);
		end;
	if TextInfo <> nil then
		with TextInfo^ do begin
				count := TextTE^^.TELength;
				err := fswrite(f, count, TextTE^^.hText^);
				err := fsclose(f);
				Convolve('TempKernel', SystemRefNum);
				err := fsdelete('TempKernel', SystemRefNum);
			end;
end;


function NewPicWindowD (name: str255): boolean;
	const
		WidthID = 5;
		HeightID = 6;
		TitleID = 8;
	var
		mylog: DialogPtr;
		item: integer;
		SaveWidth, SaveHeight: integer;
		okay: boolean;
begin
	if not macro and not OptionKeyWasDown then begin
			InitCursor;
			SaveWidth := NewPicWidth;
			SaveHeight := NewPicHeight;
			mylog := GetNewDialog(190, nil, pointer(-1));
			SetDNum(MyLog, WidthID, NewPicWidth);
			SelIText(MyLog, WidthID, 0, 32767);
			SetDNum(MyLog, HeightID, NewPicHeight);
			SetDString(MyLog, TitleID, name);
			OutlineButton(MyLog, ok, 16);
			repeat
				ModalDialog(nil, item);
				if item = WidthID then begin
						NewPicWidth := GetDNum(MyLog, WidthID);
						if (NewPicWidth < 0) or (NewPicWidth > MaxPicSize) then begin
								NewPicWidth := SaveWidth;
								SetDNum(MyLog, WidthID, NewPicWidth);
							end;
					end;
				if item = HeightID then begin
						NewPicHeight := GetDNum(MyLog, HeightID);
						if (NewPicHeight < 0) or (NewPicHeight > MaxPicSize) then begin
								NewPicHeight := SaveHeight;
								SetDNum(MyLog, HeightID, NewPicHeight);
							end;
					end;
			until (item = ok) or (item = cancel);
			if item = ok then
				name := GetDString(MyLog, TitleID);
			DisposDialog(mylog);
			if NewPicWidth < 32 then
				NewPicWidth := 32;
			if odd(NewPicWidth) then
				NewPicWidth := NewPicWidth + 1;
			if NewPicHeight < 16 then
				NewPicHeight := 16;
			if item = cancel then begin
					NewPicWidth := SaveWidth;
					NewPicHeight := SaveHeight;
					exit(NewPicWindowD);
				end;
		end; {if not macro}
	NewPicWindowD := NewPicWindow(name, NewPicWidth, NewPicHeight);
end;


procedure PlotSurface;
	var
		hend, vend, h, v, DataWidth, DataHeight, i: integer;
		htemp, vtemp, ivalue: integer;
		skip, DataLeft, DataRight, DataTop, DataBottom: integer;
		hLoc, vLoc, hMin, hMax, vMin, vMax, MinIValue, MaxIValue: integer;
		hstart, vstart, dh, dv, hbase, vbase, vscale, nPlotLines, CalValue: extended;
		peak, MaxPeak, hinc, vinc, nLines, MinCValue, MaxCValue: extended;
		poly: PolyHandle;
		SaveInfo, PlotInfo: InfoPtr;
		aLine: LineType;
		MaskRect: rect;
		AutoSelectAll, ApplyLUT: boolean;
		table: LookupTable;
		StartTicks: LongInt;

	procedure FindVinc;
	begin
		with PlotInfo^.PicRect do begin
				vstart := 5.0 + MaxPeak - dv * DataWidth;
				skip := round(DataHeight / ((bottom - vstart - 5.0) / vinc));
				if skip = 0 then
					skip := 1;
				nPlotLines := DataHeight / skip;
				vinc := (bottom - vstart - 5.0) / nPlotLines;
				vinc := vinc / 0.95;
				repeat
					vinc := vinc * 0.95;
					hinc := vinc / 2.0;
				until (5.0 + hinc * nPlotLines + dh * DataWidth) < right;
			end;
	end;

begin
	if NotRectangular or NotInBounds then
		exit(PlotSurface);
	StopDigitizing;
	DisableDensitySlice;
	SetForegroundColor(BlackIndex);
	SetBackgroundColor(WhiteIndex);
	SaveInfo := Info;
	if not NewPicWindowD('Surface Plot') then begin
			KillRoi;
			exit(PlotSurface)
		end;
	PlotInfo := info;
	info := SaveInfo;
	AutoSelectAll := not Info^.RoiShowing;
	ShowWatch;
	if AutoSelectAll then
		SelectAll(true);
	if TooWide then
		exit(PlotSurface);
	with info^ do
		ApplyLUT := ((LutMode = GrayScale) or (LutMode = CustomGrayscale)) and (not IdentityFunction);
	if ApplyLUT then
		GetLookupTable(table);
	Measure;
	UndoLastMeasurement(true);
	with results do begin
			MinIValue := MinIndex;
			MaxIValue := MaxIndex;
		end;
	if ApplyLut then begin
			MinIvalue := table[MinIValue];
			MaxIvalue := table[MaxIValue];
		end;
	MinCValue := 10e100;
	MaxCValue := -10e100;
	for i := MinIValue to MaxIValue do begin
			ivalue := i;
			if ApplyLUT then
				ivalue := table[ivalue];
			calValue := cvalue[i];
			if calValue < minCValue then
				minCValue := calValue;
			if calValue > maxCValue then
				maxCValue := calValue;
		end;
	WhatToUndo := NothingToUndo;
	with results do
		if (MaxValue - MinValue) <> 0.0 then
			vscale := (255.0 / (MaxValue - MinValue)) * 0.5
		else
			vscale := 0.5;
	with info^.RoiRect do begin
			DataLeft := left;
			DataRight := right;
			DataTop := top;
			DataBottom := bottom;
			DataWidth := DataRight - DataLeft;
			DataHeight := DataBottom - DataTop;
		end;
	dh := (0.65 * PlotInfo^.PicRect.right) / DataWidth;
	dv := -0.4 * dh;
	hstart := 5.0;
	vinc := 2.0;
	MaxPeak := (MaxCValue - MinCValue) * vscale * 0.5;
	FindVinc; {First estimate}
	MaxPeak := MaxPeak * 2.0;
	hmin := DataRight + round(MaxPeak / dv);
	if hmin < 0 then
		hmin := 0;
	vmax := DataTop + round(MaxPeak / vinc);
	if vmax > DataBottom then
		vmax := DataBottom;
	MaxPeak := 0.0;
	vloc := DataTop;
	skip := 3;
	repeat
		hloc := hmin;
		repeat
			ivalue := MyGetPixel(hloc, vloc);
			if ApplyLUT then
				ivalue := table[ivalue];
			calValue := cvalue[ivalue];
			peak := (calValue - MinCValue) * vscale + (DataRight - hloc) * dv - (vloc - DataTop) * vinc;
			if peak > MaxPeak then
				MaxPeak := peak;
			hloc := hloc + skip;
		until hloc > DataRight;
		vloc := vloc + skip;
	until vloc > vmax;
	FindVinc;
	v := DataTop;
	StartTicks := TickCount;
	SetPort(GrafPtr(PlotInfo^.osPort));
	PenNormal;
	repeat
		hmax := 0;
		vmin := 9999;
		poly := OpenPoly;
		hbase := hstart;
		vbase := vstart;
		Info := SaveInfo;
		GetLine(DataLeft, v, DataWidth, aLine);
		info := PlotInfo;
		if ApplyLUT then
			ApplyTableToLine(@aLine, table, DataWidth);
		MoveTo(round(hbase), round(vbase - vscale * (cvalue[aLine[0]] - MinCValue)));
		for i := 0 to DataWidth - 1 do begin
				hbase := hbase + dh;
				vbase := vbase + dv;
				hLoc := round(hbase);
				vLoc := round(vbase - vscale * (cvalue[aLine[i]] - MinCValue));
				LineTo(hloc, vloc);
				if hloc > hmax then
					hmax := hloc;
				if vloc < vmin then
					vmin := vloc;
			end;
		LineTo(round(hbase), round(vbase));
		LineTo(round(hstart), round(vstart));
		LineTo(round(hstart), round(vstart - vscale * (cvalue[aLine[0]] - MinCValue)));
		hmin := round(hstart);
		vmax := round(vstart);
		ClosePoly;
		ErasePoly(poly);
		FramePoly(poly);
		KillPoly(poly);
		SetRect(MaskRect, hmin, vmin, hmax, vmax);
		UpdateScreen(MaskRect);
		hstart := hstart + hinc;
		vstart := vstart + vinc;
		v := v + skip;
	until (v >= DataBottom) or CommandPeriod;
	ShowTime(StartTicks, SaveInfo^.RoiRect, '');
	if CommandPeriod then
		beep;
	info^.changes := true;
end;


procedure MakeSkeleton;
{This table-driven parallel thinning routine is based on an algorithm}
{by Zhang and Suen(CACM, March 1984, 236-239). There is}
{an entry in the table for each of the 256 possible 3x3 neighborhood}
{configurations. An entry of '1' means delete pixel on first pass, '2' means}
{delete pixel on second pass, and '3' means delete on either pass. There is a}
{routine in 'user.p' that will draw all 256 neighborhoods.}
	const
		s999 = '01234567890123456789012345678901';
		s000 = '00030033003130330000000030203033';
		s032 = '00000000300000003000000030003022';
		s064 = '00000000000000000000000000000000';
		s096 = '30000000200020003000000030003020';
		s128 = '03330013000000010000000000000001';
		s160 = '31000000000000002000000000000000';
		s192 = '33130013000000010000000000000000';
		s224 = '3301000100000000330100002200200';
	var
		table: FateTable;
		s: str255;
		i, pass: integer;
begin
	s := concat(s000, s032, s064, s096, s128, s160, s192, s224);
	for i := 0 to 254 do
		table[i] := ord(s[i + 1]) - ord('0');
	table[255] := 0;
	pass := 0;
	repeat
		PixelsRemoved := 0;
		filter(skeletonize, pass, table);
		pass := pass + 1;
		if not CommandPeriod then
			filter(skeletonize, pass, table);
		pass := pass + 1;
	until (PixelsRemoved = 0) or CommandPeriod;
end;


procedure DoErosion;
	var
		i: integer;
		t: FateTable;
begin
	for i := 0 to BinaryIterations - 1 do begin
			filter(Erosion, i, t);
			if CommandPeriod then
				leave;
		end;
end;


procedure DoDilation;
	var
		i: integer;
		t: FateTable;
begin
	for i := 0 to BinaryIterations - 1 do begin
			filter(Dilation, i, t);
			if CommandPeriod then
				leave;
		end;
end;


procedure DoOpening;
	var
		i: integer;
		t: FateTable;
begin
	for i := 0 to BinaryIterations - 1 do begin
			filter(Erosion, i, t);
			if CommandPeriod then
				exit(DoOpening);
		end;
	for i := 0 to BinaryIterations - 1 do begin
			filter(Dilation, i + BinaryIterations, t);
			if CommandPeriod then
				exit(DoOpening);
		end;
end;

procedure DoClosing;
	var
		i: integer;
		t: FateTable;
begin
	for i := 0 to BinaryIterations - 1 do begin
			filter(Dilation, i, t);
			if CommandPeriod then
				exit(DoClosing);
		end;
	for i := 0 to BinaryIterations - 1 do begin
			filter(Erosion, i + BinaryIterations, t);
			if CommandPeriod then
				exit(DoClosing);
		end;
end;

procedure SetBinaryCount;
	var
		TempCount: integer;
		Canceled: boolean;
begin
	TempCount := GetInt('Neighborhood Pixel Count(1-8):', BinaryCount, Canceled);
	if Canceled then
		exit(SetBinaryCount);
	if (TempCount >= 1) and (TempCount <= 8) then begin
			BinaryCount := TempCount;
			BinaryThreshold := BinaryCount * 255
		end
	else
		beep;
end;

procedure SetIterations;
	var
		TempIterations: integer;
		Canceled: boolean;
begin
	TempIterations := GetInt('Number of Iterations:', BinaryIterations, Canceled);
	if Canceled then
		exit(SetIterations);
	if (TempIterations >= 1) and (TempIterations < 100) then
		BinaryIterations := TempIterations
	else
		beep;
end;


procedure ChangeValues (v1, v2, v3: integer);
      {Changes all pixels in the current selection with a value in the range v1 to v2 to a value of v3.}
	var
		i, value: integer;
		table: LookupTable;
begin
	for i := 0 to 255 do begin
			value := i;
			if (value >= v1) and (value <= v2) then
				value := v3;
			table[i] := value;
		end;
	ApplyTable(table);
end;


procedure DoPropagate (MenuItem: integer);
      {Copies the current Look-Up Table, spatial calibration, or density calibration to all open windows.}
	var
		TempInfo: InfoPtr;
		i: integer;

	procedure CopyLUTInfo;
	begin
		with info^ do begin
				TempInfo^.RedLUT := RedLUT;
				TempInfo^.GreenLUT := GreenLUT;
				TempInfo^.BlueLUT := BlueLUT;
				TempInfo^.ColorStart := ColorStart;
				TempInfo^.ColorEnd := ColorEnd;
				TempInfo^.nColors := nColors;
				TempInfo^.LutMode := LUTMode;
				TempInfo^.cTable := cTable;
				TempInfo^.FillColor1 := FillColor1;
				TempInfo^.FillColor2 := FillColor2;
				TempInfo^.FillColor1 := FillColor1;
				TempInfo^.SaveFill1 := SaveFill1;
				TempInfo^.SaveFill2 := SaveFill2;
			end;
	end;

	procedure CopySpatialCalibration;
		var
			SaveInfo: InfoPtr;
	begin
		with info^ do begin
				TempInfo^.xSpatialScale := xSpatialScale;
				TempInfo^.ySpatialScale := ySpatialScale;
				TempInfo^.PixelAspectRatio := PixelAspectRatio;
				TempInfo^.RawspatialScale := RawspatialScale;
				TempInfo^.ScaleMagnification := ScaleMagnification;
				TempInfo^.Units := Units;
				TempInfo^.UnitsID := UnitsID;
				TempInfo^.FullUnits := FullUnits;
				TempInfo^.changes := true;
				TempInfo^.SpatiallyCalibrated := SpatiallyCalibrated;
			end;
		SaveInfo := Info;
		Info := TempInfo;
		UpdateTitleBar;
		Info := SaveInfo;
	end;

	procedure CopyDensityCalibration;
		var
			SaveInfo: InfoPtr;
	begin
		with info^ do begin
				TempInfo^.DensityCalibrated := DensityCalibrated;
				TempInfo^.ZeroClip := ZeroClip;
				TempInfo^.fit := fit;
				TempInfo^.nCoefficients := nCoefficients;
				TempInfo^.Coefficient := Coefficient;
				TempInfo^.UnitOfMeasure := UnitOfMeasure;
				TempInfo^.changes := true;
			end;
		SaveInfo := Info;
		Info := TempInfo;
		UpdateTitleBar;
		Info := SaveInfo;
	end;

begin
	for i := 1 to nPics do begin
			TempInfo := pointer(WindowPeek(PicWindow[i])^.RefCon);
			case MenuItem of
				1: 
					CopyLUTInfo;
				2: 
					CopySpatialCalibration;
				3: 
					CopyDensityCalibration;
			end; {case}
		end;
	WhatToUndo := NothingToUndo;
end;

procedure DoArithmetic (MenuItem: integer; constant: extended);
	var
		table: LookupTable;
		i: integer;
		tmp: LongInt;
		LogScale: extended;
		Canceled: boolean;
begin
	canceled := false;
	if not macro then
		case menuItem of
			AddItem: 
				constant := GetReal('Constant to add:', 25, Canceled);
			SubtractItem: 
				constant := GetReal('Constant to subtract:', 25, Canceled);
			MultiplyItem:  begin
					constant := GetReal('Constant to multiply by:', 1.25, Canceled);
					if constant < 0.0 then begin
							PutMessage('Constant must be positive.');
							exit(DoArithmetic);
						end;
				end;
			DivideItem:  begin
					constant := GetReal('Constant to divide by:', 1.25, Canceled);
					if constant <= 0.0 then begin
							PutMessage('Constant must be nonzero and positive.');
							exit(DoArithmetic);
						end;
				end;
			LogItem:  begin
					constant := 0.0;
					LogScale := 255.0 / ln(255.0);
				end;
		end; {case}
	if Canceled then
		exit(DoArithmetic);
	for i := 0 to 255 do begin
			case MenuItem of
				AddItem: 
					tmp := round(i + constant);
				SubtractItem: 
					tmp := round(i - constant);
				MultiplyItem: 
					tmp := round(i * constant);
				DivideItem: 
					tmp := round(i / constant);
				LogItem: 
					if i = 0 then
						tmp := 0
					else
						tmp := round(ln(i) * LogScale);
			end;
			if tmp < 0 then
				tmp := 0;
			if tmp > 255 then
				tmp := 255;
			table[i] := tmp;
		end;
	ApplyTable(table);
end;


procedure AutoThreshold;
       {Iterative thresholding technique, described originally by Ridler & Calvard in}
       {"PIcture Thresholding Using an Iterative Selection Method", IEEE transactions}
       { on Systems, Man and Cybernetics, August, 1978. }
	var
		AutoSelectAll, SaveRedirectFlag: boolean;
		index, MovingIndex, level: integer;
		tempSum1, tempSum2, tempSum3, tempSum4, result: extended;
begin
	AutoSelectAll := not info^.RoiShowing;
	if AutoSelectAll then
		SelectAll(false);
	SaveRedirectFlag := RedirectSampling;
	RedirectSampling := false;
	if info^.RoiType = RectRoi then
		GetRectHistogram
	else
		GetHistogram;
	RedirectSampling := SaveRedirectFlag;
	OptionKeyWasDown := OptionKeyDown;
	if not OptionKeyWasDown then begin
              {Default is to set to these to null so erased areas won't be included in the threshold }
			Histogram[0] := 0;
			Histogram[255] := 0;
		end;
	with Results do begin {From ComputeResults}
			MinIndex := 0;
			while (histogram[MinIndex] = 0) and (MinIndex < 255) do
				MinIndex := MinIndex + 1;
			MaxIndex := 255;
			while (histogram[MaxIndex] = 0) and (MaxIndex > 0) do
				MaxIndex := MaxIndex - 1;
			if (MinIndex >= MaxIndex) then begin
					level := 128;
					ShowMessage(concat('Threshold=', Long2Str(level)));
					EnableThresholding(level);
					exit(AutoThreshold);
				end;
			MovingIndex := MinIndex;
			repeat
				tempSum1 := 0;
				tempSum2 := 0;
				tempSum3 := 0;
				tempSum4 := 0;
				for index := MinIndex to MovingIndex do begin
						tempSum1 := tempSum1 + index * Histogram[index];
						tempSum2 := tempSum2 + Histogram[index];
					end;
				for index := (MovingIndex + 1) to MaxIndex do begin
						tempSum3 := tempSum3 + index * Histogram[index];
						tempSum4 := tempSum4 + Histogram[index];
					end;
				Result := (tempSum1 / TempSum2 / 2) + (tempSum3 / tempSum4 / 2);
				MovingIndex := MovingIndex + 1;
			until ((MovingIndex + 1) > result) or (MovingIndex > (MaxIndex - 1));
			level := Round(result);
			EnableThresholding(level);
			ShowMessage(concat('Threshold=', Long2Str(level)));
		end; {with}
end;


procedure AutoDensitySlice;
	var
		AutoSelectAll: boolean;
		sigmak1k2, sigmax, nsum: real;
		i, j, maxk1, maxk2, temp: integer;
		musubt, omegak1, omegak2, muk1, muk2: real;
		part1, part2, part3: real;
		intermed1, intermed2, intermed3: real;
begin
	ResetGrayMap;
	AutoSelectAll := not info^.RoiShowing;
	if AutoSelectAll then
		SelectAll(false);
	if info^.RoiType = RectRoi then
		GetRectHistogram
	else
		GetHistogram;
	maxk1 := 0;
	maxk2 := 0;
	musubt := 0.0;
	nsum := 0.0;
	for i := 1 to 254 do begin
			nsum := nsum + histogram[i];
		end;
	for i := 1 to 254 do begin
			musubt := musubt + (i * (histogram[i] / nsum));
		end;
	sigmak1k2 := 0.0;
	sigmax := 0.0;
	omegak1 := 0.0;
	muk1 := 0.0;
	for i := 1 to 253 do begin
			temp := i + 1;
			omegak2 := 0.0;
			muk2 := 0.0;
			omegak1 := omegak1 + (histogram[i] / nsum);
			muk1 := muk1 + (i * (histogram[i] / nsum));
			if omegak1 > 0.0 then begin
					for j := temp to 254 do begin
							omegak2 := omegak2 + (histogram[j] / nsum);
							muk2 := muk2 + (j * (histogram[j] / nsum));
							if omegak1 * omegak2 * (1.0 - omegak1 - omegak2) > 0.0 then begin
									part1 := ((omegak1 * muk2) - (omegak2 * muk1)) * ((omegak1 * muk2) - (omegak2 * muk1));
									intermed1 := omegak2 * omegak1;
									part2 := ((omegak1 * (musubt - muk2)) - (muk1 * (1 - omegak2))) * ((omegak1 * (musubt - muk2)) - (muk1 * (1 - omegak2)));
									intermed2 := omegak1 * (1 - omegak1 - omegak2);
									part3 := ((omegak2 * (musubt - muk1)) - (muk2 * (1 - omegak1))) * ((omegak2 * (musubt - muk1)) - (muk2 * (1 - omegak1)));
									intermed3 := omegak2 * (1 - omegak1 - omegak2);
									if intermed1 * intermed2 * intermed3 > 0.0 then begin
											sigmak1k2 := part1 / intermed1 + part2 / intermed2 + part3 / intermed3;
										end;
								end;
							if sigmak1k2 > sigmax then begin
									maxk1 := i;
									maxk2 := j;
									sigmax := sigmak1k2;
								end;
						end;
				end;
		end;
	SliceStart := maxk1;
	SliceEnd := maxk2;
end;


procedure FixColors;
    {Because Image always sets LUT entries 0 and 255 to white and black respectively we need to map}
    {pixels with values of 0 or 255 to the nearest matching color in the other 254  LUT entries.}
	var
		i, index2, match0, match255: integer;
		table: LookupTable;

	procedure BestMatch (index1: integer; var match: integer);
		var
			i: integer;
			rdiff, gdiff, bdiff: LongInt;
			diff, mindiff: extended;
	begin
		match := index1;
		mindiff := 10e10;
		if index1 = 0 then
			index2 := 1
		else
			index2 := 254;
		with info^ do
			for i := 1 to 254 do begin
					rdiff := bsr(cTable[index1].rgb.red, 8) - bsr(cTable[index2].rgb.red, 8);
					gdiff := bsr(cTable[index1].rgb.green, 8) - bsr(cTable[index2].rgb.green, 8);
					bdiff := bsr(cTable[index1].rgb.blue, 8) - bsr(cTable[index2].rgb.blue, 8);
					diff := sqrt(sqr(rdiff) + sqr(gdiff) + sqr(bdiff));
					if diff < mindiff then begin
							match := index2;
							mindiff := diff;
						end;
					if index1 = 0 then
						index2 := index2 + 1
					else
						index2 := index2 - 1;
				end;
	end;

begin
	BestMatch(0, match0);
	BestMatch(255, match255);
	table[0] := match0;
	for i := 1 to 254 do
		table[i] := i;
	table[255] := match255;
	ApplyTable(table);
end;


procedure GetDItemRect (d: DialogPtr; item: integer; var r: rect);
	var
		iType: integer;
		ignore: handle;
begin
	GetDItem(d, item, itype, ignore, r)
end;


procedure DrawPopUpText (str: str255; r: rect);
begin
	TextFont(SystemFont);
	if (str = '+') or (str = 'Ð') or (str = 'Ö') then begin
			TextSize(24);
			MoveTo(r.left + 13, r.bottom - 2);
		end
	else begin
			TextSize(12);
			MoveTo(r.left + 13, r.bottom - 5);
		end;
	DrawString(str);
end;


procedure ImageMathUProc (d: DialogPtr; item: integer);
     {User proc for Image Math dialog box}
	var
		str: str255;
		VersInfo: str255;
		r: rect;
begin
	SetPort(d);
	GetDItemRect(d, item, r);
	DrawDropBox(r);
	case item of
		OpItem:  begin
				GetItem(ImageMathOpsMenuH, ord(CurrentMathOp) + 1, str);
				DrawPopUpText(str, r);
			end;
	end;
end;


procedure SetUProc (d: DialogPtr; item: integer; pptr: handle);
	var
		itype: integer;
		r: rect;
		h: handle;
begin
	GetDItem(d, item, itype, h, r);
	SetDItem(d, item, itype, pptr, r);
end;


procedure DoImageMath;
	const
		ScaleItem = 10;
		OffsetItem = 11;
		ResultItem = 12;
	var
		d: DialogPtr;
		item, i, MenuItem: integer;
		r: rect;
		str: str255;
begin
	InitCursor;
	d := GetNewDialog(200, nil, pointer(-1));
	SetUProc(d, Src1Item, @ImageMathUProc);
	SetUProc(d, Src2Item, @ImageMathUProc);
	SetUProc(d, OpItem, @ImageMathUProc);
	repeat
		if item = OpItem then begin
				setport(d);
				GetDItemRect(d, item, r);
				MenuItem := PopUpMenu(ImageMathOpsMenuH, r.left, r.top, ord(CurrentMathOp) + 1);
				case MenuItem of
					1: 
						CurrentMathOp := AddMath;
					2: 
						CurrentMathOp := SubMath;
					3: 
						CurrentMathOp := MulMath;
					4: 
						CurrentMathOp := DivMath;
					5: 
						CurrentMathOp := AndMath;
					6: 
						CurrentMathOp := OrMath;
					7: 
						CurrentMathOp := XorMath;
					8: 
						CurrentMathOp := MaxMath;
					9: 
						CurrentMathOp := MinMath;
					10: 
						CurrentMathOp := CopyMath;
				end;
				DrawDropBox(r);
				GetItem(ImageMathOpsMenuH, ord(CurrentMathOp) + 1, str);
				DrawPopUpText(str, r);
			end;
		ModalDialog(nil, item);
	until (item = ok) or (item = cancel);
	DisposDialog(d);
	if item = cancel then
		exit(DoImageMath);
end;


end.